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^S^ • The hydrodynamic limit for a kinetic model of chemotaxis is investigated. The limit 

C^ i equation is a non local conservation law, for which finite time blow-up occurs, giving rise 

to measure- valued solutions and discontinuous velocities. An adaptation of the notion of 

duality solutions, introduced for linear equations with discontinuous coefficients, leads to 

^\j . an existence result. Uniqueness is obtained through a precise definition of the nonlinear 

^ I flux as well as the complete dynamics of aggregates, i.e. combinations of Dirac masses. 

)^ • Finally a particle method is used to build an adapted numerical scheme. 
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^ 1 Introduction 

Kinetic frameworks have been investigated to describe the chemotactic movement of cells in the 
presence of a chemical substance since in the 80 's experimental observations showed that the 
motion of bacteria (e.g. Escherichia Coli) is due to the alternation of 'runs and tumbles'. The 
so-called Othmer-Dunbar-Alt model [H [I2l [201 122] describes the evolution of the distribution 
function of cells at time t, position x and velocity v, assumed to have a constant modulus c > 0, 
as well as the concentration S{t, x) of the involved chemical. A general formulation for this 
model can be written as 

dtfe + V-VJe = -I {T[Se]{v' ^ v) ^v') - T[Se]{v ^ v')Uv)) dv' , 



-ASe + Se = peit,x):= / fe{t,X,v)dv. 

'\v\=c 



The second equation describes the dynamics of the chemical agent which diffuses in the domain. 
It is produced by the cells themselves with a rate proportional to the density of cells p and 
disappears with a rate proportional to 5*. The transport operator on the left-hand side of the first 
equation stands for the unbiased movement of cells ('runs'), while the right-hand side governs 
'tumbles', that is chemotactic orientation, or taxis, through the turning kernel T[S]{v' — )■ f), 
which is the rate of cells changing their velocity from w' to u. 

The parameter e corresponds to the time interval of information sampling for the bacteria, 
usually e ^ 1, and when it goes to zero, one expects to recover the collective behaviour of the 
population, that is a macroscopic equation for the density p(t, x) of cells. Such derivations have 
been proposed by several authors. When the taxis is small compared to the unbiased movement 
of cells, the scaling must be of diffusive type, so that the limit equations are of diffusion or 
drift- diffusion type, see for instance [9J for a rigorous proof. In [T^ [21], the authors show that 
the classical Patlak-Keller-Segel model can be obtained in a diffusive limit for a given smooth 
chemoattractant concentration. 

We focus here on the opposite case, that is when taxis dominates the unbiased movements. 
This is accounted for in the model by the choice of the scaling in equation fll.ip . Moreover, we 
consider positive chemotaxis, which means that the involved chemical is attracting cells, and 
therefore is called chemoattractant. The model has been proposed in [lOj, several works have 
been devoted to the mathematical study of this kinetic system. Existence of solutions has been 
obtained for various assumptions on the turning kernel in [HI EJ [HI [TB]. Numerical simulations 
of this system are proposed in [27]. The limit problem is usually of hyperbolic type, see for 
instance [T3| [23| [2^ for a hyperbolic limit model which consists in a conservation equation for 
the cell density and a momentum balance equation. 

It is not difficult to obtain the following formal hydrodynamic limit to equation fll.ip . more 
precisely on the total density of particles p = lim^ p^: 

dtp + dw,{a[S]p) = 0, -A,,5 + S = p. (1.2) 

Here the macroscopic velocity a[S] depends on the chemoattractant concentration S through 
the turning kernel. This system of equations has been obtained in [10], with a rigorous proof in 
the two-dimensional setting for a fixed smooth S, and therefore a bounded density p. The aim 
of this paper is to obtain rigorously this limit for the whole coupled system. Severe difficulties 
arise then mainly due to the lack of estimates for the solutions to the kinetic model when e goes 
to zero and consequently to the very weak regularity of the solutions to the limit problem. 

It turns out that the limit equation is in some sense a weakly nonlinear conservation equation 
on the density p. Indeed the expected velocity field depends on p, but through S, and therefore 
in a non local way. Actually it can be written as a variant of the so-called aggregation equation, 
for which blow-up in finite time is evidenced (see e.g. [3]), leading to measure- valued solutions. 
In this respect, this equation behaves also like linear equations with discontinuous coefficients. 
In particular Dirac masses can arise, this is the mathematical formulation of the aggregation of 
bacteria. Therefore S is no longer smooth, and a major difficulty in this study will be to define 
properly the velocity field a = a[S'] and the product ap. 

The viewpoint of the aggregation equation has been extensively studied by Carrillo et al. j8] 
through optimal transport techniques. Existence and uniqueness are obtained in a very weak 
sense, and the dynamics of aggregates is also given. We propose here another approach, based 
on the notion of duality solutions, as introduced in the linear case by Bouchut and James |1]. 



The main drawback is that presently we have to restrict ourselves to the one-dimensional case, 
since the theory in higher dimensions is not complete yet (see [6j). The approach proposed by 
Poupaud and Rascle [26], which coincides with duality in the 1-d case, could also be explored. 
Notice however, that the properties of the expected velocity field a in the two-dimensional case 
are not obvious either. 

More precisely, we propose to proceed in a similar way as in [5j, where the nonlinear system 
of zero pressure gas dynamics is interpreted as a system of two linear conservation equations 
coupled through the definition of the product. This last point turns out to be crucial in order 
to obtain a proper uniqueness result for the system ( II. 2p . In this work, the product ap will be 
defined thanks to the limiting fiux of the kinetic system (11. ip (see also [16] for another application 
of the same idea). As we shall see, this is closely related to the dynamics of aggregates, that is 
combinations of Dirac masses, which refiect some kind of collective behaviour of the population. 
Finally, an important application of this aggregate dynamics is the development of a numerical 
scheme, based on a particle method. The motion and collapsing of Dirac masses is clearly 
evidenced. 

The paper is organized as follows. In Section [5] we precisely state the model. Section [2] is 
devoted to the notion of duality solutions, and contains the main results of this article. Some 
technical properties which will be useful for the rest of the paper are given in Section |1[ Then 
we investigate in Section [5| the proof of the existence and uniqueness result of duality solution 
for system (I2.8p - (l2.10p stated in Theorem 13.91 In Section [H] we prove the rigorous derivation of 
the hydrodynamical system from the kinetic system. Finally, the dynamics of aggregates and 
the numerical scheme for the limit equation are described in the last section, where numerical 
illustrations are also provided. 

2 Modelling 

From now on we focus on the one dimensional version of the problem, so that x G M. We first 
recall the main assumptions leading to the kinetic equation, next we proceed to the formal limit. 

2.1 Kinetic model 

In this work, cells are supposed to be large enough to sense the gradient of the chemoattractant 
instantly. Therefore the turning kernel takes the form (independent on v) 

T[S]{v' ^v) = ^{v'd^S). (2.1) 

The function $ is the turning rate, obviously it has to be positive. More precisely, for attractive 
chemotaxis, the turning rate is smaller if cells swim in a favourable direction, that is v-V xS > 0. 
Thus $ should be a non increasing function. A simplified model for this phenomenon is the 
following choice for $: we fix a positive parameter a, a mean turning rate 0o > and take 

$(x) = 0o(l + 0(x)), (2.2) 

where (t> is an odd function such that 

IT T* <r — (y 

gC°°(M), 0'<O, 0(x) = <; '7 .^ ' (2.3) 

if X > a, 



where < A < 1 is a given constant. 

Now since the transport occurs in M the set of velocities is ^ = {— c, c}, and the expression 
of the turning kernel simplifies in such a way that (11.11) rewrites 

dtfe + vdj, = ^m-vd,S)M-v) - ^vd^S)fe{v)), V E V. (2.4) 

- d,,S, + S, = pe = fe{c) + fe{-c). (2.5) 

The existence of weak solutions in a L^ setting for a slightly different system in a more 
general framework has been obtained for instance in [71 [15]. Concerning precisely this model, 
we refer to [27] for the existence theory in any space dimension. Notice that no uniform L°° 
bounds can be expected. The reader is referred to [27] for some numerical evidences of this 
phenomenon, which is the mathematical translation of the concentration of bacteria. This is 
some kind of "blow-up in infinite time", which for e = leads to actual blow-up in finite time, 
and creation of Dirac masses. Moreover the balanced distribution vanishing the right hand side 
of (12. 4 p depends on Ss] thus the techniques developed e.g. in [9] cannot be applied. 

2.2 Formal hydro dynamic limit 

We formally let e go to assuming that S'^ and /^ admit a Hilbert expansion 

fe = fo + ^fi^ , Se = So + eSi^ 

Multiplying (12. 4p by e and taking e = 0, we find 

<l>{-cd,So)fo{-c) = <^{cd,So)fo{c). (2.6) 

Summing equations (12. 4p for c and — c, we obtain 

dtifeic) + U-c)) + cd^ifeic) - /,(-c)) = 0. (2.7) 

Moreover, from equation (12. 6p we deduce that 

/o c - /o -c = — — ... /o c + /o -c . 
$(-09^.50) + ^[cd^So) 

The density at equilibrium is defined by p := /o(c) + fo{—c). Taking £ = in (12. 7p we finally 
obtain 

dtp + d^{a{dr,So)p) =0, 

where a is defined by 

a{d^So) = c— — . . . = -c(f){cd^So), 

and we have used (12. 2p for the last identity. Notice that a is actually a macroscopic quantity, 
since it is the simplified formulation of 

^(^ o. IvvHvd.So)dv 

Jy^{vdxSo)dv 



in the one-dimensional context. 

We couple this equation with the limit of the elliptic problem (I2.5P for the chemoattractant 
concentration, so that, in summary, and dropping the index 0, the formal hydrodynamic limit 
is the following system 

dtp + dM{d^S)p) = ^, (2.8) 

a{d,S) = -c(P{cd^S), (2.9) 

-d^^S + S = p, (2.10) 
complemented with the boundary conditions 

p(t = 0,x) = p™*(x), lim p(t, x) = 0, lim 5(t,a;) =0. (2.11) 

x— >-±oo a:— >-±oo 

We now give the precise formulation of the limit system in terms of aggregate equation. 
Noticing that a solution to (l2.1Up has the explicit expression 

S{t,x) = K*p{t,.){x), where ir(x) = ^e-l^', (2.12) 

the macroscopic conservation equation for p (12. 8p can be rewritten 

dtp + d^{a{d^,K * p)p) = 0. 

When a is the identity function, this is exactly the so-called aggregation equation, and since the 
potential is non-smooth, blow-up in finite time is expected. We refer the reader to e.g. [3, 8j, 
and [T^ in the context of chemotaxis. 

Similar problems were encountered for instance in [18j, where the authors investigate the 
high field limit of the Vlasov-Poisson-Fokker-Planck model in one space dimension. The limit 
system is a scalar conservation law coupled to the Poisson equation, and a proper definition of 
the product is needed to pass to the limit. This definition has been extended in two dimensions 
by Poupaud [2S] using defect measures but losing uniqueness. 



3 Duality solutions 

3.1 Notations 

Let Cq{Y, Z) be the set of continuous functions from Y to Z that vanish at infinity and Cc{Y, Z) 
the set of continuous functions with compact support from Y to Z. All along the paper, we 
denote A4ioc{^) the space of local Borel measures on R. For p G Aiioc we denote by |p|(M) its 
total variation. We will denote A^b(M) the space of measures in Aiioc(^) whose total variation 
is finite. From now on, the space of measure-valued function A^ft(]R) is always endowed with 
the weak topology a{Mb,Co). We denote Sm ■= C{[0,T];Mb(^) - a{Mb,Co)). 

We recall that if a sequence of measure {pn)nm in AitiM) satisfies sup„gjij |/Xn|(R) < +oo, 
then we can extract a subsequence that converges for the weak topology a{Aib, Co). 

The coupled system (I2.8|) - (l2.9p - (l2.10j) is interpreted in this context as a linear conservation 
equation (12. 8p . the velocity b of which depends on the solution S to the elliptic equation (I2.10p . 



b = a{dxS). This actually means that equation (12. 8p is somehow nonlinear. One convenient 
tool to handle such conservation equations 

dtP + dxipp) = 0, h being a given function, (3.1) 

whose solutions eventually are measures in space, is the notion of duality solutions, introduced 
in [4J. 

3.2 Linear conservation equations 

Duality solutions are defined as weak solutions, the test functions being Lipschitz solutions to 
the backward linear transport equation 

dtp + hit, x)d,p = 0, p{T, .)=p^ e Lip(M). (3.2) 

A key point to ensure existence of smooth solutions to (13. 2p is that the velocity field has to be 
compressive, in the following sense. 

Definition 3.1 We say that the function h satisfies the so-called one-sided Lipschitz condition 
(OSL condition) if 

dxb{t, .) < (3{t) for [i G L-'^(0,T) in the distributional sense. (3-3) 

A formal computation shows that dt{pp) + dx[b{t,x)pp] = 0, and thus 

d 



^^» , p{t,x)p{t,dx)] =0, (3.4) 

which defines the duality solutions for suitable p's. It is now quite classical that (13. 3p ensures 
existence for (13. 2p . but not uniqueness, which is of great importance here to obtain stability 
results and make a convenient use of (13. 4p . 

Therefore, the corner stone in the construction of duality solutions is the introduction of the 
notion of reversible solutions to (13. 2p . A complete statement of the definitions and properties 
of reversible solutions would be too long in the present context, so that merely a few hints are 
given. Let C denote the set of Lipschitz continuous solutions to (13. 2p . and define the set of 
exceptional solutions: 

S = <p (z C such that p = 0>. 
The possible loss of uniqueness corresponds to the case where £ is not reduced to zero. 

Definition 3.2 We say that p E C is a reversible solution to ^3.^) if p is locally constant on 
the set 

Ve = |(t,x) e [0,T] X M; 3 pe G ^, Peit.x) ^ O}. 

This definition leads quite directly to the uniqueness results of |1]. It turns out that the 
class of reversible solutions is also stable by perturbations of the coefficient h. 

We now restrict ourselves to those p's in (13. 4p . More precisely, we state the following defini- 
tion. 



Definition 3.3 We say that p e Sm '■= C{[0,T];Mb{M.) - a{Mb,Co)) is a duality solution 
to Ii3.1\) if for any < r < T, and any reversible solution p to liS.^) with compact support in 

X, the function t i— )■ / p{t,x)p{t,dx) is constant on [0,r]. 

Jr 

Remark 3.4 A similar notion of duality solution for the transport equation is available dtU + 
hdxU = 0, and p is a duality solution of Ii3.1\) iff u = J p is a duality solution to transport 
equation (see ^). 

We shall need the following facts concerning duality solutions. 
Theorem 3.5 (Bouchut, James ^) 

1. Given p° G A^{,(]R), under the assumptions (|2I3J, there exists a unique p G Sm, duality 
solution to /i3. 1\) . such that p(0, .) = p° ■ 

Moreover, if p° is nonnegative, then p{t, ■) is nonnegative for a.e. t>0. And we have the 
mass conservation 

|p(t,-)|(K) = |p°|(M), for a.e. tG]0,T[. 

2. Backward flow and push-forward: the duality solution satisfies 

VtG [O,T],V0gCo(M), [ (p{x)p{t,dx)= f ^{X{t,0,x))p%dx), (3.5) 

Jr Jr 

where the backward flow X is defined as the unique reversible solution to 
dtX + b{t, x)dxX = zn]0,s[xR, X{s,s,x) = x. 

3. For any duality solution p, we define the generalized flux corresponding to p by h^p = 
—dtU, where u = J pdx. 

There exists a bounded Borel function b, called universal representative of b, such that 
b = a almost everywhere, and for any duality solution p, 

dtp + dx{bp) = in the distributional sense. 

4. Let (bn) be a bounded sequence in L°°{]0,T[xM.), such that bn^b in L°°(]0,T[x]R) — w-k. 
Assume dxbn < an(^)j where {an) is bounded in L^{]0,T[), dxb < a E L^{]0,T[). Consider 
a sequence {pn) G Sm of duality solutions to 

dtPn + dxibnPn) = m ]0,T[xM, 

such that Pn{0, .) is bounded in Aih(R), and p„(0, .) ^^ p° E Aib{W). 
Then pn ^ p in Sm, where p G Sm is the duality solution to 

dtp + dx{bp) = m ]0,r[xM, p(0, .) = p°. 
Moreover, bnPn -^ bp weakly in A^5(]0,T[x]R). 



The set of duality solutions is clearly a vector space, but it has to be noted that a duality 
solution is not a priori defined as a solution in the sense of distributions. However, assuming 
that the coefficient b is piecewise continuous, we have the following equivalence result: 

Theorem 3.6 Let us assume that in addition to the OSL condition fl3.3p . b is piecewise con- 
tinuous on ]0,T[x]R where the set of discontinuity is locally finite. Then there exists a function 
b which coincides with b on the set of continuity of b. 

With this b, p E Sm is a duality solution to (13. ip if and only if dtp + dx{bp) = in V{M.). 
Then the generalized flux b^p = bp. In particular, b is a universal representative ofb. 

This result comes from the uniqueness of solutions to the Cauchy problem for both kinds of 
solutions (see Theorem 4.3.7 of |1]). 

3.3 Main results 

We are now in position to give the definition of duality solutions for the limit system (I2.8p - (l2.10p . 

Definition 3.7 We say that (p, 5) G C([0,T]; A^b(R)) x C{[0,T];W^'°^) is a duality solution 
to d^ - d^TTU]) if there exists b E L°°((0,r) x M) and a E L]^^{Q,T) satisfying dj < a in V, 
such that 

1. for alio <ti<t2<T 

dtp + dxibp) = in the sense of duality on ]ti, t2[, 

2. ( 12. 9p is satisfied in the weak sense: 

VVgC^(M), VtG [0,T], I {dxSdx^ + S^){t,x)dx= f ip{x) p{t,dx), 

Jr J 

3. b = a{dxS) a.e. 

Remark 3.8 For S m C([0,T]; 1^^'°°) and 4> as in (Q, we have ai^d^S) E C([0,T]; L°°(M)). 
Therefore equation (12. 8 P is meaningful in the duality sense. The key property is then the one- 
sided Lipschitz condition. 

Unfortunately, Definition 13.71 does not ensure uniqueness, as we shall evidence in Section 
O This is due to the fact that the product a{dxS)p is not properly defined yet. Indeed the 
relevant definition of this product relies on a proper definition of the fiux of the system, which 
we introduce now. Let A be an antiderivative of a such that v4(0) = 0, we set 

J = -dM{dxS)) + a{dxS)S. (3.6) 

This choice is justified first since this definition holds true when S is regular. Indeed we have 
dx{A{dxS)) = a{dxS)dxxS, so that we can write J = a{dxS){—dxxS + S) = a{dxS)p. On the 
other hand, a more physical reason relies on the fact that the above J is the correct fiux for the 
kinetic model, and passes to the limit when e goes to zero, see Section [61 
We can now establish the following uniqueness theorem: 



Theorem 3.9 Let us assume that p*™ > is given in A^b(M). Then, for allT > there exists 
a unique duality solution (p, S) with p > of (12. 8p -( I2.10p which satisfies in the distributional 
sense: 

dtp + dj = 0, (3.7) 

where J is defined in (13. 6p . It means that the universal representative in Theorem \3.5\ satisfies 

bp = J, in the sense of measures. 

Moreover, we have p = X^p™* where X is the backward flow corresponding to a{dxS). 

The second resuh concerns the rigorous proof of hydrodynamical hniit for the kinetic model. 
Let (/g, Se) be a solution of the system (I2.4p - (l2.5p . complemented with null boundary condition 
at infinity and with the following initial data: 

/,(0,-,-) = /r, (3.8) 

such that p*"* = Ve * P*™ where t]^ is a mollifier and p*"* is given in A^f,(]R). We recall that 
for fixed e > 0, there exists (/g, Se) such that /^ belongs to C([0,T] x R x y) and therefore 

SeeC{[0,T];C^{R)), see [7], or [27] in the present context. 

Theorem 3.10 Let us assume that p*™ > is given in A^;,(M). Let (/g, S^) be a solution to the 
kinetic-elliptic equation (I2.4p - (l2.5p with initial data (13. 8p . Then, as £ — )■ 0, (/jjS'e) converges 
in the following sense: 

p, := /,(c) + fe{-c) - p in Sm ■■= C{[0, T]; A^,(M) - a{Mt, Co)), 

S,^S m C{[0,T];W^''^{R))-weak, 

where (p, S) is the unique duality solution of the system (12. 8p - (I2.10p satisfying 

bp = J, in the sense of measures. 

4 Properties of S 

We gather in this section a set of properties for the solution S to (I2.10p that will be used 
throughout the paper. 

4.1 One-sided estimates 

The estimates presented in this part rely only on equation (I2.10p . 

Lemma 4.1 Let p G C([0,T], A^{,(M)). Then the solution S of equation ()2.10p satisfies 

1. p>0^5>0 

2. one-sided estimate: d^xS < S if and only if p > 

3. for all p e [l,+oo], S e C{[0,T],LP{R)) and d^S e C{[0,T],Lp{R)) 



9 



Proof. The first two items are easy consequences of tlie expression (I2.12p for tlie first one, of 
tlie equation (I2.10p for tlie second. For tlie third item, from convolution properties, we have for 

any p G [1, +00] 

||5(t, .)||lp{r) = T^lle"''' *p{t,-)\\Lp{R) < \p{t,.)\{R)-\\e-\-^\\LP{R) = - sup |p(t, ■)|(M), 
^ z / te[o,T] 

where |p|(M) stands for the total mass of the nonnegative measure p. We proceed in the same 
way for d^S. □ 

As mentioned above, the key point to use the duality solutions is that the velocity field 
satisfies the OSL condition (I3.3p . 

Lemma 4.2 Let p G Sm- Then the coefficient a{dxS) defined by fl2.9p - fl2.10p satisfies the OSL 
condition fl3.3p if and only if p > 

Proof. Straightforward computations lead to 

dMd.S)) = -c^cP'{cd,S)d^^S. 

With (12.1 op and since is a nonincreasing function, we deduce from the one-sided estimate of 
Lemma 14.11 

5.(o(5.5))<max{cl0'|Uoo^,O}. 

We conclude thanks to the bound on S in L°°. q 

Finally, we turn to a convergence result for a sequence of such functions S. 

Lemma 4.3 Let {pn)nen be a sequence of measures that converges weakly towards p in Sm o-s 
n goes to +00. Let Sn{t,x) = {K * Pn{t, ■)){x) and S{t,x) = {K * p{t, ■)){x), where K is defined 
in fl2.12p . Then when n — )■ +00 we have 

d^Snit.x) — > d^S{t,x) /or a.e. t G [0,T], xeM, 
dxSn(t,x) -^ dxS{t,x) in Lf^weak — *. 

Proof. The proof of this result is obtained by regularization of the convolution kernel (see 
Lemma 3.1 of [17j). ^ 



4.2 Entropy estimates 

In this subsection, we consider now that (p, S) satisfy fl3.7p -f l3l6|) in the sense of distributions. 
We prove first that S satisfies a nonlinear nonlocal equation. Next, following the strategy of 
|19] . we prove that the above one-sided estimate implies some kind of entropy inequality for 
dxS. 

Lemma 4.4 Assume {p,S) G C([0,T]; A^fc(M))xC([0,r]; H/i'°°) satisfy (^1\ -^M. thend^S G 
C([0,T],Li(M)) nL~([0,T],5K(M)) and S is a weak solution of 

dtS - dxK * dMidxS)) + dxK * {a{dxS)S) = 0. (4.1) 
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Proof. We have p E Sm and d^^S = S - p. Then d^S E C([0,T], ^^(M)) nL°°([0,T], 5y(M)). 
We recall that we have S = K * p where K{x) = |e~'^''. Thus taking the convolution by K of 
f l3.7p -f l3r6|) . we get that 5* is a weak solution of (14. ip . q 

Lemma 4.5 Let S be a weak solution in C{[0,T];W^'^(R)) of (14.1 p with initial data S**™. 
We assume moreover that d^S belongs to L°°{[0,T]; BV(R)) and that the one-sided estimate 
dxxS < S holds in the distributional sense. Then for any twice continuously difjerentiable 
convex function t] we have 

dtvidxS) + dMd.S)) - r]'{dxS)a{dxS)S + r]'{dxS)[K * {-d^Aid^S) + a{dxS)S)] < 0, (4.2) 

where the entropy flux q is defined by 

PX 

q(x) = / Ti\y)a{y) dy. 



Proof. From Lemma r4.4[ S satisfies (14. ip . By differentiation, and using the property dxxK = 
K — 5q, we get 

dtdxS + dxA{dxS) - a{dxS)S + K * {-dxA{dxS) + a{dxS)S) = O. (4.3) 

Consider a sequence of mollifiers C,n{x) = n(^{nx), with n E N, ( E C^(R), C ^ aiid 
/k C(3^) dx = 1. We set Sn = Cn * S. Then we have 

dtdxSn + dx{A{dxS) * Cn) - Cn * {a{dxS)S) + K*Cn* {-dxA{dxS) + a{dxS)S) = 0. 

We define the commutators i?„ and Qn as follows: 

A{dxS)*Cn = A{dxSn)+Rn{t,x), 

Qn{t,x) = -Cn*{a{dxS)S) + K*Cn*{-dxA{dxS) + a{dxS)S) 

+ a{dxSn)Sn - K * {-dxA{dxSn) + a{dxSn)Sn), 

SO that the regularized solution satisfies 

dtdxSn + dx{A{dxSn) + Rn) ' a{dxSn)Sn + K * (-9,A(9,5„) + a{dxSn)Sn) + Q„ = 0. (4.4) 

Let us consider 77 a twice continuously differentiable convex function and let q be the corre- 
sponding entropy flux. Multiplying equation (14.40 by ri'{dxSn), we get 

dtvidxSn) + dx{q{dxSn) + v'{dxSn)Rn) + Hn = -v'{dxSn)Qn + Rndx{ri'{dxSn)), (4.5) 

where 

Hn ■■= -v'idxSn)aidxSn)S^ + v'idxSn)[K * {-dxA{dxSn) + aidxS^)Sn)]. 

Due to properties of the convolution product, we have 

Rn ^0, Qn^O in ^L((0> 00) xR), l<p< +00, 
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so that in the sense of distribution, we have straightforwardly 

d.W{d,Sn)Rn) ^ 0, v\d,Sn)Qn ^ 

and 

H^^H:= ~ri'id,S)a{d,S)S + r]'id,S)[K*i-d,Aid,S) + aid^S)S)], 

which is precisely the desired term in the limit equation. Now we deal with the term Rndx{ri'{dxSn)) 
on the right-hand side, and we notice that Rn > thanks to the Jensen inequality and the con- 
vexity of A. Therefore, since rj is convex, we have 

RndAr]'{d,Sn)) = Rnr]"{d,Sn)d,,Sn < i?„r/"(9,5„)5„, 

where we have used the one-sided estimate d^xSn < Sn to obtain the last inequality. Since S'„, 
is bounded in L°° independently of n, we can pass to the limit in this last identity thanks to 
the Lebesgue dominated convergence theorem to get 

Rnri"{dA)Sn -> in LL((0, oo) x M). 

Finally, letting n going to +00 in fl4.5p . we deduce that (14.21) holds in the distributional sense. 
D 

Remark 4.6 This equation relies strongly on the definition of the flux J in (13. 6p . This fact 
has already been noticed by the authors in 116], which can be viewed as a particular case of 
the one studied in this paper by replacing the elliptic equation f l2.10p for S by the Poisson 
equation —d^xS = p. In this case, the product of a{dxS) by p is naturally defined by a{dxS)p = 
—dxA{dxS), so that equation on S corresponding to (14. 3p is given by 

dtdxS + dxA{dxS) = 0. 

This equation is a nonlinear hyperbolic conservation law which is local, contrary to (14. Sp . There- 
fore uniqueness is ensured by entropy conditions. Since dxS is monotonous (—dxxS = p > 0), 
this can be formulated as a chord condition on A (see [S]). If in addition A is convex or concave 
(i.e. if a is non- decreasing or non-increasing), this selects only increasing or decreasing shocks. 

5 Existence and uniqueness for the hydrodynamical prob- 
lem 



In this Section, we focus on the proof of Theorem 13. 9[ which can be split in 3 steps. The 
first one consists in obtaining the dynamics of aggregates, or in other words of combinations of 
Dirac masses. Next we obtain the existence of duality solutions in the sense of Definition 13.71 
by proving first that aggregates define such a solution, then proceeding to the general case by 
approximation. This is exactly the same strategy as for the pressureless gases in [5j. Finally, 
uniqueness follows from a careful definition of the flux of the equation. In this respect, we first 
underline with an example that Definition 13.71 as it stands does not give uniqueness, and how 
the proper definition of the flux singles out a unique solution. 
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Indeed, let us consider (I2.8p -( 12?T0|) with boundary condition ( I2.1ip where the initial datum 
is assumed to be a Dirac mass in 0: p*™ = 6q. We have that {60, K * 60) is a solution to 
f l2.8p -f l2TT0|) with initial data 6q. Actually, the pair 

p,{t,x) = 5,,(^t){x); S,{t,x) = K*p,{t,x) = ie-l^-^^WI. (5.1) 

turns out to define a solution in the sense of duality in Definition 13.71 for several choices of 
curves Xi with Xi(0) = 0. Set bi{t,x) = a{dxSi){t,x), and notice first that, according to 
Remark 13. 4^ pi is a duality solution if ui := J pidx = H{x — Xi{t)) is a duality solution of 
the transport equation. Now, from Lemma [4.21 61 satisfies the OSL condition, therefore Ui is a 
duality solution of the transport equation as soon as it is solution in the sense of distributions. 
As detailed in [4j, Section 3, this holds true only if u satisfies some admissibility conditions, 
namely, the characteristics of the velocity field have to enter the discontinuity on both side. 
Since lim^_^^+ 6i(x) = a(— 1/2) and lim^_j.^- bi{x) = a(l/2), the velocity of the shock should 
satisfy a(l/2) > x'i{t) > a(— 1/2), which furnishes an infinity of solution. 

For any of the previous solutions, the generalized fiux given by Theorem 13.51 -3 is biApi = 
—dtUi = —x[{t)6xj^{t)- On the other hand, let us compute the flux J defined by (13.60 . For 
simplicity, we set here a = in the definition (12.30 of $. With this convention, we get 

f -Ac, X < xi(t), 1 f -Ace^"''^^*^ x<xi(t), 

a{d.Si){t,x) = \ A{d.Si){t,x) = -\ 

[ Ac, x>xi{t), 2 1^ _Ace"-'+''^^*'', x > xi{t). 

Obviously we have J = 0, so that the condition ap = J selects x[{t) = 0, which finally implies 
Xi = since a;i(0) = 0. 

5.1 Dynamics of aggregates 

Let us first consider the motion of aggregates. We assume that p^""' is given by a finite sum of 
Dirac masses: p™* = ^27=1 ''^i^x° where xj < a;2 < ■ ■ ■ < a;° and the rrii-s are nonnegative. We 
look for a couple (p„, S'„) solving in the distributional sense dtPn + d^Jn = where the fiux J„ 
is given by (13. 6p and Sn solves (I2.10p . We recall that it means that Sn = K * pn where K is 
defined in (I2.12p . Let us set pn{t,x) = Y17=i^^i^Xi(t)- Such a function is a solution in the sense 
of distributions of (13. 7p if the function Un defined by 

Unit.x) := j pndx = '^miH{x - Xi{t)), (5.2) 

where H denotes the Heaviside function, is a distributional solution to 

dtUn - d,A{d,Sn) + a{d.,Sn)Sn = 0. (5.3) 

We have 



2 



n 



d,Sn{t,x) = -^!^ Sign (x-x,(t))e-l^-^''WI. (5.4) 
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Straightforward computations prove that we have in the distributional sense 

n 

d,A{d,Sn) = a{d,Sn)Sr, + J^iMd.Sn)],^., (5-5) 

j=l 

where [/]a;. = f{xf) — f{x~) is the jump of the function / at Xj. Injecting f l5.2p and ( 15. Sp in 
(15. 3p . we find 

n n 

i=l j=l 

Thus the dynamics of aggregates is finally given by 

mix'iit) = -[A{d^Sn)]x,(t), for z = 1, . . . , n. 

We complement this system of ODEs by the initial data Xj(0) = x°. More precisely, recalling 
that K(x) = |e~'^', using (15. 4p this latter system can be rewritten : 

— ^ + 2J mjdxK{xj -Xi)\ - AI — -^ + 2J rnjda:K{xj - Xi) \ . (5.6) 

J9^i J \ i+i J 

Recall that, from the definition of the coefficient a in (12. 9p with ( 12. 3p . a is nondecreasing and 
odd, so that A is a convex function. This implies that for i = 1, . . . , n — 1, x^ > a^^+i, therefore, 
aggregates can collapse in finite time but an aggregate cannot split. This is a direct consequence 
of the fact that we are considering positive chemotaxis, i.e. a is nondecreasing. If there exists a 
time ti for which we have for instance Xiit^ = Xj+i(ti), then the dynamics for t > ti is defined 
as above except that we replace rrii by rrii + rrii+i and Xi{t) = Xi+i(t) for t > ti. Moreover A is 
even, then when n = 1, we have x[ = and Xi{t) = a:?. Thus if aggregates collapse such that 
they form a single aggregate of mass ^^ rrii, then this aggregate does not move for larger times. 

5.2 Existence of duality solutions 

We have constructed (p„, Sn) which is a solution of ( I3.7p -( l3l6l) -( l2.10p in the distributional sense 
for the given initial data p^™. We recall the following result due to Vol'pert [28j (see also [2j): 
if u belongs to BV{R) and / G C^(M) with /(O) = 0, then v = f ou belongs to BV(R) and 

^lu with 7„ = f'{u) a.e. such that (/ o u)' = J^u'. 

Together with the fact that A is an antiderivative of a such that A{0) = 0, this result implies 
that there exists a function a„ such that 

J„ := -d:r{A{dr,Sn)) + a{dxSn)Sn = anPn, and a„ = a{dxSn) a.e. 
Thus pn is a solution in the distributional sense of 

dtpn + dxijinPn) = 0. 

Moreover, we deduce from (15. 4p that a((9^5'„) is piecewise continuous with the discontinuity 
lines defined by x = Xi, i = 1, . . . ,n. We can apply Theorem 13.61 which gives that p„ is a 
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duality solution and that a„ is a universal representative of a{dxSn)- Then the flux is given by 

a{dxSn)APn = Jn- 

Let us yet consider the case of any initial data p*™ G Aib{W). We approximate p*™ by 
Pn"* — Y17=i''^i^x° with p™ ^^ p*"* in Aib{^)- By the same token as above, we can construct 
a solution (p„, Sn = K * pn) with p„(t = 0) = p^'^* = X^iLi ^^i^x°i which solves in the sense of 
duality 

dtPn + dx{a{dxSn)pn) = 0, 

in the sense of distributions 

dtpn + dx-Jn = 0, J„ = -(9^A((9^5'„) + a{dxSn)Sn, 
and which satisfies 

Moreover, since dxSn is bounded in L°° uniformly with respect to n by construction, we can 
extract a subsequence of {a{dxSn))n that converges in L°° — weak* towards h. Since from 
Lemma 14. 2[ a{dxSn) satisfies the OSL condition, we deduce from Theorem 13.51 4) that, up 
to an extraction, p„ ^ p in Sm and anPn -^ Q-P weakly in A^;,(]0,T[x]R), p being a duality 
solution of the scalar conservation law with coefficient h. With Lemma 14. 3[ we deduce that 
dxSn -^ dxS a.e., it implies in particular that Jn ^ J '■= —dxA{dxS) + a{dxS)S in V{W) and 
that a{dxSn) — > a{dxS) a.e. By uniqueness of the weak limit, we have b = a{dxS). Moreover 
J = ap a.e. and p satisfies then (13. 7p . Then (p, S) is a solution as in Theorem 13. 9 [ this concludes 
the proof of the existence. 

5.3 Uniqueness of solutions 

Let us consider yet the study of the uniqueness. As shown above. Definition 13.71 is not sufficient 
to ensure uniqueness. Therefore, we will use the fact that we have a duality solution p that 
satisfies (13. 7p in T>'{[0, T] x M) with the initial data p™* and with the fiux J given by (13.61) . This 
equation leads to the non-local evolution equation on S (14. ip as stated in Lemma 14.41 

Another key point is the one-sided estimate dxxS < S. In fact, if we consider for instance 
pim _ g^ then it is obvious that p = is a solution of (l2.8p -( l2rT0|) . However, if we allow p to 
be nonpositive, i.e. if the corresponding chemoattractant concentration S does not satisfy the 
one-sided estimate dxxS < S, then we can build a simple example of non-uniqueness. Indeed 
we have that 

p{t,x) = 6-xi{t){x) - 26o{x) + (5xi(t)(x) 

is a duality solution of (I2.8p - (l2.10p which satisfies (13. 7p . provided a;i(0) = and (15. 6p is satisfied. 
This readily gives 

x[{t) = A(l + e-- + le-2-) -A{-^ + e-- + le'^-). 
Here by convexity of A, we have x[ > 0. 
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Theorem 5.1 Let Si and S2 be two weak solutions in C{[0,T]]W^'^(M)) of (14. ip with ini- 
tial data 5*^"' and 5*2"* respectively. If we assume moreover that d^Si and 8x82 belongs to 
L°°([0,T]; i?y(]R)) and that the one-sided estimate 



holds in the distributional sense. Then there exists a nonnegative constant C such that 

Mo Oil ^-^ /^ 1 1 CT'tTii cT'i'ni \ \ 

ll'Jl — '-'2||l°°([0,T];W^i.1(]R)) — ^ 11*^1 "'-'2 IIVKi'i(K)- 

Proof. We start from the entropy inequahty (14. 2p of Lemma H75l Using standard regularization 
arguments, it is weU-known that we can apply this inequality to the family of Kruzkov entropies 
Vniu) = \u — k\. Then, the doubling of variables technique developed by Kruzkov allows to 
justify the following computation. Assume 5*1 and 5*2 are two weak solutions of (14. ip . then in 
the distributional sense, we have 

dtmSi -82)1+ dxi signidxSi - dxS2){A{dxSi) - Aid^))) < 
sign(9,^i - 8x82) {dxK * {A{dxSi) - A{dxS2)) + aiSi - 02^2 - K * {a^Si - 02^2)), 

where we denote ai = a{dxSi) and 02 = 0(9^.5*2). Integrating with respect to x and using the 
properties of the convolution product, we deduce 

^ f \dx{Si - ^2)1 dx < \\dxK\\^ [ \A{dxSi) - A{dxS2)\dx + (1 + \\K\\^) [ {a^S^ - a2S2\ dx. 

The function a being regular, we have 

^ / \dx{Si - S2)\dx <Co [ mS, - S2)\dx + Ci /" 1^1 - S2\ dx. (5.7) 

"* Jr Jr Jr 

In the same way as for equation (14. ip . this leads to 
d 



^^ , ,5i - ^2! dx<C2 I mSi - 52)1 dx + Cs I \Si - S2\ dx. (5.8) 

Summing (15.80 and (15. 7p . we deduce that there exists a nonnegative constant C such that 

— IPi — 5*21114/1.1(18) < C\\Si — 52||vKi'i(R)- 

Applying the Gronwall Lemma allows to conclude the proof. q 

Proof of uniqueness in Theorem 13.91 Let us assume that we have two duality solutions 
(pi. Si) and (p2, ^2) such as in Theorem 13.91 Therefore, from Lemma SiU Si and ^2 are weak so- 
lutions of (14. ip . Using Theorem 15. H we conclude that 5i = ^2. Thus pi = K * Si = K * S2 = P2- 

n 
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6 Convergence for the kinetic model 

In this section we investigate the rigorous derivation of ( 12.8p -( l2TT0l) from the microscopic model 
(12. 4p . First we state some estimates on the moments of the solution of the kinetic problem. 

Lemma 6.1 Let (f^jSg) be a solution of the kinetic problem fl2.4p - fl2.5p . Then for all t G [0,T] 
and all 6 > we have 

[ [ \v\^fedxdv= |t;|'=|p"*|(M), k en. 
Jr Jv 

Proof. Since v eV = {—c,c}, \v\ is constant therefore 

/ / \v\''fsdxdv= |f|^ / psdx. 
Jr Jv Jr 

The result follows then directly from the mass conservation in (12. 4p . q 

Proof of Theorem [3A0l Let (fe, Ss) be a solution of i^^-I^B)- For fixed e > 0, we have 
/, e C([0,T] xRxV). Define p, := fyf^dv, J, := JyVf.dv and a{d^S,) = -c0(c9,5,). We 
can rewrite the kinetic equation (12. 4p as 

dtfe + vdj, = ^m-vd,S,)p, - 2/,). 
Taking the zeroth and first order moments, we get 

5ip, + 9,J, = 0, (6.1) 

2 
dtJe + v^d^p^ = -{a{d^Se)pe - Je)- (6.2) 

From ([61]), we deduce that Vt G [0,r], \p,{t, ■)|(M) = |p™|(R). Therefore, for all t G [0,T] 
the sequence (pe(t, ■))£ is relatively compact in A1b(]R) — cr(A^{,(]R), Co(M)). Moreover, there 
exists Ue G L°°([0, T], _BV(R)) such that p^ = dxUe- From (16. ip . we get that dtUe = —Je and 
thanks to Lemma l6.ll we deduce that m^ is bounded in Lip([0,T], L^(R)). This implies the 
equicontinuity in t of {ps)e- Thus the sequence {ps)e is relatively compact in Sm and we can 
extract a subsequence still denoted (pe)e that converges towards p in Sm- 

We recall that Se{t,x) = [K * pe{t,-)){x) where K{x) = |e~'^L Denoting S{t,x) := {K * 
p{t,-)){x), since p G S^, we have d^S G L°°{[0,T]; BV{Wj). From Lemma [4.3[ the sequence 
{dxSs)e converges in L°°w — * and a.e. to dxS as e goes to 0. Lemma [4.21 ensures that both 
aldxSe) and a{dxS) satisfy the OSL condition. 

From (I6.ip -( l6l2|) . we have in the distributional sense 

dtps + dxia{dxSe)ps) = dxia{dxSe)pe - Je) = ^dxidtJe + v'^dxpe) = Re- (6.3) 

Now, for all tp G C^((0,T) x M), we deduce from Lemma [6.11 



{dtJe + v'^dxPe)dxipdxdt 



< \v\\pmR)\\^t^x^P\\L^ + \v\'\pmRm 
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This implies that the hmit in the distributional sense of the right-hand side Re of fl6.3p vanishes. 
Now we multiply equation (12.51) by a^d^S;;) and use again the antiderivative A of a to obtain 

a{d,Se)p, = -dMid-,Se)) + a{d^Se)Se, (6.4) 

so that we can rewrite the conservation equation (16. 3p as follows, in ©'(M): 

dtp, + d, {-d,A{d,Se) + a{d,S,)Se) = ^-dMJe + v^d.pe). (6.5) 

Taking the limit e — )■ of equation (16. 5p in the sense of distributions, we get 

dtp + d, {-d,A{d,S) + aid,S)S) = in V'{R), (6.6) 

where S(t,x) = {K * p{t,-)){x). Therefore the pair {p, S) satisfies (I3.7p - (l3.6p . In addition, p 
is nonnegative as a limit of nonnegative measures, so that Lemma 14.11 implies the one-sided 
estimate d^xS < S. Thus we are in position to apply Lemma [4.41 and Theorem 15. ![ which give 
uniqueness for S, and consequently for p. Therefore the whole sequence p^ converges to p in 
<Sm- We recall that we have chosen the initial data such that p^"* = Ve * P*"* where 77^ is a 
mollifier. Therefore p^™ ^ p™^ in MbW - (r{Mb{^),Co{R)). 

Thus we have constructed a solution that satisfies (16. 6 P in the distributional sense, in other 
words, we have defined a solution of the problem (I2.8p - (l2.10p thanks to its flux. A natural 
question is to know whether we can define a velocity corresponding to this flux. From the 
theory of duality solutions (see Theorem 13. 5p . it boils down to show that the above constructed 
solution is a duality solution. From Vol'pert calculus [2S] we infer the existence of as such that 
as = a{dxS) a.e. and 

dxiAidxS)) = asdxxS. 

Therefore 

— dx{A{dxS)) + a{dxS)S = asP a.e. , with as = a^d^S) a.e. (6.7) 

Using equation (16. 6p we have in the distributional sense 

dtp + dx{asp)=0. (6.8) 

However, we have proved in Section 15.31 that such a solution is unique. We deduce that the 
solution (p, 5*) obtained by the hydrodynamical limit above is the duality solution of Theorem 
13.91 It concludes the proof of Theorem 13.101 ^ 

Remark 6.2 In the proof above, the macroscopic flux J defined in (13. 6p appears to be the limit 
of the microscopic flux J^. Indeed from (16. 2 p and (16. 4p we deduce that, in the distributional 
sense, 

J,^J:= -dxAidxS) + aidxS)S. 

This natural definition of the flux allows to get the uniqueness of the solutions of the coupled 
system (I2.8p - (l2.10p thanks to equations (I4.ip - (l4.3p . Such a technique to establish the hydrody- 
namic limit has been proposed in fT8^. But the authors do not state that their limit is a duality 
solution and do not define a velocity and therefore a flow corresponding to their flux. In the 
limit of the Vlasov-Poisson-Fokker-Planck system, this result has been investigated in JTOj - 



7 Numerical issue 

7.1 Finite time of collapse 

Before focusing on the numerical simulations, let us clarify the dynamics of the model. In the 
case of n Dirac masses, rrij > for i = 1, . . . , n, located at positions Xi < ■ ■ ■ < x^, we recall 
that the time evolution is governed by system ( 15. 6p : 

-Y + z2 ^j9xK{xj -Xi)\ - A — -^ + 22 mjd^K{xj - Xi) \ , (7.1) 

for i = 1, . . . ,n, where we recall that A is an antiderivative of a such that A{0) = 0. We deduce 



that for all t > 0, and for i = 1, . . . ,n, 



— ^ + X] ^j9xK{xj -Xi),Y + Yl ^j9xK{xj - Xi) 1 
such that x'i{t) = a{'~fi{t)). 
Proposition 7.1 Let us assume that there exists n G N* such that 

n 

r{x) = J2m%o{x), 

i=l 

with m° > 0, for i = 1, . . . ,n. We assume in addition that a is a nondecreasing and odd real 
function. Then the duality solution p of Theorem \3.9\ has the following properties : 

1. Ifn = 1, xi{t) = X? for all t > 0. Then p{t) = p^™ for all t > 0. 

2. For i = 1, . . . , n — 1 , x'iit) > x'^^iit) therefore Xj+i — Xj < x^^i — x°. 

3. There exists c* G [x'j',x^] and T* > such that p{t,x) = 6c*{x) for all t > T* . 

Proof. The first point is a direct consequence of the even character of A whereas the second 
point comes from the convexity of A. Let us then prove the third point. By convexity of the 
function A and with (17. ip . we have 

m^x[>Al— + Y,^e^' "0 -^I-^ + Et'^ ) ^°' 

and 

As for (17. 2p . we can rewrite these last inequalities as : 

x[{t) > a(7i(0)) > 0, x'^< a(7„(0)) < 0. 

We deduce that there exists a time T* > such that all masses collapse for t = T* in a single 
Dirac mass. ^ 
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Remark 7.2 Notice that we have in addition the following estimate forT*: 

T*<(4-x?)/(a(7i(0))-a(7„(0))). 

Corollary 7.3 Let us assume that < p*™ G Cc(M) with compact support [0,L]. Let us denote 
p the duality solution of Theorem \3.9[ with initial data p*"*. Then there exists c* G [0, L] and 
T* >0 such that p{t,x) = 5c*{x) for all t > T* . 

Proof. Let us approximate p*™ by 

n 



i=l 

IxO 



with x° = {i — l)L/n, for i = l,...,n and m^ = J o'^^ p*"*((ix). From Proposition 17. ![ we 

i 

deduce that there exists c* G [0, L] and T* > such that the duahty solution of Theorem 
13.91 with initial data p"" is such that pn{t,x) = 5c* for all t > T*. Moreover, we have T* < 
L/(a(7"(0)) — a(7^(0))) where we recall that 

n n 

and 

n 71 

TTl Tfl 

- Y^ _^(k0-n)Lln ^ ^n(Q) < ^0 _ ^ ^gO-n)L/n_ (74) 

By stability results on duality solutions in Theorem 13.51 (see also subsection 15. ip . we deduce 
that p„ ^ p in Sm as n — )■ +00. Taking the limit in (17. 3p and (17. 4p . we deduce by continuity 
of p*"* that 

L 



lim 7i"(0)= / p™(s)e-"da; 



and 

lim 7^(0) = - / p*™(a;)e-^+"da;. 

Moreover, since p*™ is continuous with compact support in [0, L] we have p*™(0) = p*™(L) = 0. 
We deduce that the sequence (T^)„gN* is bounded. Thus there exists a time T* independent of 
n such that Pn(t) = 5c* for all t > T*. Taking the limit when n — )■ +00, we conclude that there 
exists c & [0) -^] such that p{t) = 6c for all t > T*. q 

Remark 7.4 Taking a = Id, therefore A{x) = x'^/2, we deduce from (17. ip that 

x'i = ^rajd^K^Xj - Xi). 

We recover the dynamics of the aggregation equation as noticed by Carrillo et al. in J^. These 
authors prove in particular the concentration in finite time of the total mass in the center of 
mass. In the framework of the present work, which is focused on applications to chemotaxis, 
a is not assumed to be the identity function, so that the center of mass is not conserved. A 
numerical evidence of this phenomenon will be proposed in the last subsection of this paper. 
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7.2 Discretization 

The numerical resolution of system fl2.8p -f l2?T0|) is far from obvious. A first naive idea consists 
in applying a standard splitting method where we treat separately the scalar conservation law 
f l2.8p and the elliptic equation fl2.1Up . It turns out that such a scheme is unable to recover the 
correct definition of the flux and therefore of the product a{dxS) by p. In particular, it leads to 
stationary Dirac masses. 

A second idea consists in solving the distributional conservation law ( 13. 7p by a finite volume 
method. It involves a discretization of the flux J on the interface of each cell of the mesh, and 
thus one could expect a correct computation of the flux, and therefore a convenient interpretation 
of the product. However, this definition of the flux involves the calculation of two derivatives 
of S. Using a centered scheme to discretize this quantity induces spurious oscillations as it 
is usually noticed for centered scheme on scalar conservation laws. We can then upwind the 
scheme depending on the sign of a{dxS) computed at previous iteration. But in doing so, we 
actually specify a value for a{dxS) in the definition of the product a{dxS) with p, and this can 
lead to capture wrong solutions. 

Next, one can think of solving the equation (14. ip on S, motivated by the fact that it plays a 
key part in the uniqueness, and that p can be recovered readily from 5*. However the equation 
is non local and its numerical resolution appears to be quite complicated and with a high 
computational cost (even in the one dimensional setting). 

Thus we prefer to use a method based on the dynamics of aggregates, detailed in Section [5?T1 
We use the principle of a particle method in which we approximate the density by a sum of Dirac 
masses. Then the motion of these pseudo-particles is approximated by discretizing system (15.60 
with an explicit Euler scheme. More precisely, let us assume that we have an approximation of 
p at time t„ = nAt, given by 

p-{x) = J2<^?i^), (7-5) 

where m^ > is the mass allocated to the pseudo-particle at the position y^ with y^ < y"^ < 
■ ■ ■ < y^n for /" G N*. Then an approximation of the potential at time t" is given by 

7" 

Using an explicit Euler scheme, we compute the new position 







+ E 



Next, we test if some pseudo-particles have collided during the time step At. If y^^^ < y"^^ for 
j > 1, then the pseudo-particles j and j + 1 have collapsed and form a unique pseudo-particle 
which has the mass m^-|-m^_,_p In this case, we decide to set this pseudo-particle at the position 
|(?/?+i + y7^^) and set m'^^'^ = rn^ + ^^?+i, moreover we have therefore 1"+^ = /" — !. Finally, 
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for given initial sequences (yf )j=i,...,7-o and {m^)i=i...jo of size J", we can construct (y") and (rn") 
of size /" as above. 

Using well-known result on the convergence of Euler scheme, we deduce that, for given initial 
data {yi)i=i,...jo, (m°)j=i jo and /°, y" defined above converges to the solution Xi(t) of fl5.6p 
when At tends to such that t„ — )■ t. Using the convergence result in Section 15.21 we deduce 
that the function p^ in (17. 5p converges in Sm to the unique duality solution of Theorem 13. 9[ 
Then the method introduced above is convergent provided we discretize the initial data p*"* in 
such a way that p°(x) := X] j=i '^i ^y? (^) converges in Aih to p*"*. Moreover, we verify easily 
that we have 

^ m° = ^ ml, and J" < J°, for all n G N, 
and that the approximation p^ of p(t„) is nonnegative. 

7.3 Numerical results 

In this Section, we present numerical simulations of model (l2.8p - fl2.1Up using the algorithm 
described above. We first approximate the initial data p*"* > 0, which is assumed to be com- 
pactly supported for numerical purpose, in the following way: we introduce a discretization 
Xj = Xq+jAx of the bounded domain which includes the compact support of p*™ and we define 

Then the sequence {y'j)j is defined by the nodes (xj) for which m}^ is not zero, and /" correspond 
to the number of i G N such that m^ is not zero. We construct then the approximation of p*"* 

by 

1° 

We present in Figure [T] the dynamics of the density p and of the chemoattractant concen- 
tration S for an initial data p*"* given by the sum of two Gaussian functions, more precisely 

As expected, we first observe the formation of two Dirac masses at the position where dxS 
initially vanishes. Then, the two aggregates collapse in the center. Looking at the time evolution, 
we notice that the first step of formation of aggregates is fast compared to the time of collapse. 
In Figure [2] we display the dynamics for an initial data given by the sum of three Gaussian 
functions: 

We observe the formation of three Dirac masses that moves according to the dynamical system 
(17. ip . They collapse then in finite time. 

Finally, as we have already noticed, we evidence that the center of mass is not fixed. For 
instance. Figure [3] represents the dynamics of the density and of the potential for an initial data 
made of one big bump with one small bump: 
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Figure 1: Dynamics of the density p (top) and of the potential S (bottom) for an initial density 
given by the sum of two Gaussian. 

The square shows the time dynamics of the center of mass. We observe that the center of mass 
at the final time is not located at the same position as at the initial time. 
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